Study of the ST2 model of water close to the liquid-liquid critical point 
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We perform successive umbrella sampling grand canonical Monte Carlo computer simulations 
of the original ST2 model of water in the vicinity of the proposed liquid-liquid critical point, at 
temperatures above and below the critical temperature. Our results support the previous work 
of Y. Liu, A.Z. Panagiotopoulos and P.G. Debenedetti [J. Chem. Phys. 131, 104508 (2009)], 
who provided evidence for the existence and location of the critical point for ST2 using the Ewald 
method to evaluate the long-range forces. Our results therefore demonstrate the robustness of the 
evidence for critical behavior with respect to the treatment of the electrostatic interactions. In 
addition, we verify that the liquid is equilibrated at all densities on the Monte Carlo time scale of 
our simulations, and also that there is no indication of crystal formation during our runs. These 
findings demonstrate that the processes of liquid-state relaxation and crystal nucleation are well 
separated in time. Therefore, the bimodal shape of the density of states, and hence the critical 
point itself, is a purely liquid-state phenomenon that is distinct from the crystal-liquid transition. 



I. INTRODUCTION 

In 1992, a numerical investigation of the equation of 
state (EOS) of the ST2 model [1] in the supercooled re- 
gion suggested the possibility of a liquid-liquid (LL) criti- 
cal point in water [5] . This initial study has subsequently 
generated a large amount of numerical and experimental 
work 3 12 . In addition to the conceptual novelty of a 
one-component system with more than one liquid phase, 
the existence of the associated LL critical point can also 
rationalize many of the thermodynamic anomalies which 
characterize liquid water (e.g. the density maximum and 
compressibility minimum), and which become more pro- 
nounced in the supercooled regime. Furthermore, the ex- 
istence of two distinct liquid phases of supercooled water 
can explain the polyamorphism which characterizes the 
glassy phase [T2TI14) . Indeed, simulations suggest that 
the low density amorphous solid form of water is similar 
to the structure of the low density liquid (LDL) phase, 
while the relaxed very-high density amorphous solid is 
related to the high density liquid (HDL) [15 . 

Evidence in support of a liquid-liquid critical point in 
water, and in other liquids with tetrahedral structure, 
has increased over time. A number of classical mod- 
els for water, including the recently developed and op- 
timized TIP4P/2005 [16], exhibit a van der Waals inflec- 
tion in their EOS at low temperature T that is evidence 
of phase coexistence between two liquid states [5J [T71 [TS] • 
The occurrence of a LL transition has also been proposed 
for silica [19j . and more recently, evidence for a LL crit- 
ical point and its associated thermodynamic anomalies 
have been presented for the Stillinger- Weber model of 
silicon [20] , 

Indeed, it is notable that the most compelling evidence 



for LL critical points has been generated in silico |21| . In 
almost all cases, LL critical points are predicted to occur 
in deeply supercooled liquids, where crystallization (in 
experiments) has so far prevented direct observation of 
such phenomenon in bulk systems. Compared to experi- 
ments, LL phase transitions are more readily observed 
in numerical studies because heterogeneous nucleation 
is not a factor, and the small system size (usually less 
than one thousand molecules) decreases the probability 
of observing the appearance of a critical crystal nucleus 
in the simulation box on the time scale of typical sim- 
ulations. Computer simulations have thus allowed the 
study of the liquid EOS under deeply supercooled condi- 
tions, on time scales longer than the structural relaxation 
time of the liquid but smaller than the homogenous nu- 
cleation time. Under these conditions, equilibrium within 
the metastable basin of the liquid free energy surface can 
be achieved without interference from crystal nucleation 
processes. 

Nonetheless, evaluations of the EOS via simulations 
in the canonical ensemble, or at constant pressure, do 
not provide a way to accurately estimate the location of 
the LL critical point found in water models, or to deter- 
mine its universality class. Only recently, in 2009, Liu, et 
al. [22] reported the first numerical investigation of ST2 
water in the LL critical region, performing simulations 
in the grand canonical ensemble for different values of T 
and of the chemical potential /i, and implementing Ewald 
sums to account for the long-range contributions to the 
electrostatic interactions. In this important contribution, 
the authors provided for the first time evidence of a den- 
sity of states that is a bimodal function of the density p, 
a necessary feature of a LL critical point. Importantly, 
the authors also showed that the fluctuations of the or- 
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der parameter (a combination of density and energy) are 
consistent with the expected shape for a critical system 
in the Ising universality class. 

More recently, Limmer and Chandler [23] have ques- 
tioned the interpretation of all previously published sim- 
ulations based on the ST2 potential, arguing that "that 
behaviors others have attributed to a liquid-liquid tran- 
sition in water and related systems are in fact reflections 
of transitions between liquid and crystal." In the case of 
the recent calculations of Liu, et al. [22], Ref. [23] pro- 
poses that "the Liu et al. result is a non-equilibrium 
phenomenon, where a long molecular dynamics run initi- 
ated from their low-density amorphous phase and run at 
constant T and P will eventually equilibrate in either the 
low density crystal or (more likely) in the higher density 
metastable liquid." It is thus of paramount importance 
to independently check the findings of Liu, et al. [22] , and 
at the same time, test whether or not the LDL phase is 
truly a disordered liquid phase characterized by a well- 
defined metastable equilibrium that is distinct from the 
crystal phase. 

In this article, we conduct these tests by carrying out 
an independent evaluation of the density of states based 
on the successive umbrella sampling technique [23]. We 
implement the original ST2 model, with the reaction field 
correction for the long range electrostatic forces, rather 
than Ewald sums, to be able to strictly compare our re- 
sults with previously published data for ST2 [25] , as well 
as to test if the LL transition is robust and independent 
of the treatment of the long range interactions. As we 
show below, we find that our results are entirely consis- 
tent with those of Liu, et al. [32], as well as with earlier 
simulation data. We further find that there is no contri- 
bution to the density of states due to crystal formation, 
confirming the distinct existence of both the HDL and 
LDL phases for T less than T c , the temperature of the 
LL critical point. 



II. MODEL AND SIMULATION METHODS 

We study the original ST2 potential as defined by Rah- 
man and Stillinger [Tj. with reaction field corrections to 
approximate the long-range contributions to the electro- 
static interactions. In the ST2 potential, water is mod- 
eled as a rigid body with an oxygen atom at the cen- 
ter and four charges, two positive and two negative, lo- 
cated at the vertices of a tetrahedron. The distances 
from the oxygen to the positive and negative charges are 
0.1 and 0.08 nm, respectively. The oxygen-oxygen inter- 
action is modeled using a Lennard- Jones potential with 
(Xlj = 0.31 nm and €lj = 0.31694 kJ/mol. We trun- 
cate this Lennard- Jones interaction at 2.5<jlj, account- 
ing for the residual interactions through standard long 
range corrections, i.e. assuming the radial distribution 
function can be approximated by unity beyond the cutoff. 
The charge-charge interactions are smoothly switched off 
both at small and large distances via a tapering function, 



as in the original version of the model pQ. 

Our grand canonical Monte Carlo (MC) algorithm is 
based on roto-translational moves, insertions, and dele- 
tions, each attempted with ratios 2:1:1. Our simulation 
box is cubic with sides of length 2 nm. The displacement 
move is accomplished by a random translation in each 
direction of up to ±0.01 nm and a random rotation of up 
to ±0.2 rad around a random direction, resulting in an 
acceptance ration of about 50%. Insertion and deletion 
moves have a much smaller acceptance ratio, of the order 
of 10~ 5 . The simulations have thus been performed for 
more than 10 10 attempted insertion/deletion moves. To 
determine the dependence of the density of states on T 
we have investigated four distinct temperatures, T = 260, 
250, 245 and 240 K. Previous numerical estimates based 
on the EOS indicate T c = 247 ± 3 K (25J [25] . 

To study the phase behavior of the system we imple- 
ment successive umbrella sampling (SUS) MC simula- 
tions [23] , from which we evaluate the probability density 
P(p) for the values of p sampled by the equilibrium sys- 
tem at fixed T, p, and volume V, and in which the num- 
ber of molecules N in the system fluctuates. In the SUS 
method, the pertinent range of p to be investigated, writ- 
ten in terms of the lower and upper number of molecules 
(respectively Ni and N u ), is divided into many small 
overlapping windows of size AN. For each window, a 
separate grand canonical MC simulation monitors how 
often a state of N particles is visited. Moreover, the 
simulations are constrained using appropriate boundary 
conditions on N, such that deletions or insertions that 
would cause N to vary outside the range assigned to 
that window are rejected [27]. The density histograms 
for each window can then be combined to obtain the full 
P{p) curve, by imposing the equality of the probability 
at the overlapping boundary. In our study Ni = 200 
(corresponding to a minimum density p — 0.75 g/cm 3 ) 
and N u = 327 (corresponding to p = 1.22 g/cm 3 ). We 
have chosen AN = 2, i.e. N is only permitted to take on 
one of two adjacent integer values within each window. 

The SUS method has a number of advantages. The 
use of narrow windows in N allows an effective sampling 
of the microstates without the use of biasing functions. 
Since the windows are independent, all the simulations 
can be run in parallel, with a gain in throughput that 
scales linearly with the number of processors employed. 
In our case, approximately 130 processors are used for 
the calculations at each T, one for each window. At 
the lowest T, more than two months of simulation time 
for each window is required for good sampling. Once 
we obtain P(p) at fixed T and p, histogram reweighting 
techniques [2S] can be applied to obtain P{p) at different 
values of p. Keeping track of the coupled density-energy 
histogram during the SUS simulations also allows us to 
estimate P(p) at different T via temperature reweight- 
ing. Finally, each window provides accurate information 
on a specific density, allowing us to compare the results 
with EOS data from previous simulations in the canoni- 
cal ensemble. 
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To facilitate comparison of our data with future stud- 
ies, in the following we will report the activity as z* = 
A x A y A z exp(/3/i)/A 3 (in units of nm" 3 ), where p is the 
chemical potential, /? is {ksT)~ l , ks is the Boltzmann 
constant, A Q = yj (2TrI a kBT) / h, I a is the principal mo- 
ment of inertia in the direction a, h is the Planck con- 
stant and A is the De Broglie wavelength, z* is the quan- 
tity that enters in the MC acceptance probability in the 
insertion-deletion moves. 



III. COMPARISON WITH PREVIOUSLY 
PUBLISHED DATA 

Before discussing the behavior of P(p), we compare our 
results for the liquid EOS as obtained from our SUS simu- 
lations, with the best available published data. In partic- 
ular, we focus on the potential energy E and the pressure 
P. The pressure is evaluated from the virial using con- 
figurations sampled at each density. Fig. [T] shows results 
from Ref. [25], obtained from molecular dynamics (MD) 
simulations, compared with our MC grand-canonical SUS 
results. At all T and p, we find excellent agreement be- 
tween these two completely independent numerical meth- 
ods. 

It is worth noting the prominent minimum observed in 
E near p — 0.83 g/cm 3 , the so-called optimal network 
density [30] . At this density the system has the pos- 
sibility to be able to satisfy all possible bonds, reaching 
at low T the ideal random tetrahedral network state. At 
lower densities, gas-liquid phase separation intervenes, 
while at larger densities closer packing prevents the sys- 
tem from satisfying the angular and distance constraints 
required to form linear hydrogen bonds between all pairs 
of molecules. 



IV. LIQUID-STATE EQUILIBRIUM 

In order to establish that the LL phase transition is 
a genuine liquid-state phenomenon, we must confirm (i) 
that all our simulations are carried out over a time scale 
that is much longer than the structural relaxation time of 
the liquid; and (ii) that the time scale for crystal nucle- 
ation is much longer than the time scale for liquid-state 
relaxation. The separation of these two time scales pro- 
vides the "window" within which the equilibrium behav- 
ior of a supercooled liquid can be defined and quantified. 

The ability of modern computer simulations to estab- 
lish liquid-state equilibrium in simulations near a LL crit- 
ical point is well documented [TBHT51 I2U1 [251 [3U] • In par- 
ticular, Ref. [30] presents the dynamical behavior of the 
same ST2 model as is studied here, as determined from 
MD simulations. In Fig. 2(a) of Ref. [30], it is shown 
that the self diffusion constant D for oxygen atoms is 
greater than 10 -8 cm 2 /s at T = 240 K for all densities 
from 0.87 to 1.2 g/cm 3 . This density range spans the 
same range within which we find bimodal behavior for 



P(p) at T = 240 K (see below). For the ST2 system, 
D > 10 -8 cm 2 /s corresponds to a range of a-relaxation 
times T a < 20 ns [3T] . Time scales well in excess of 20 ns 
are readily accessible in current MD simulations, espe- 
cially for a small system of a few hundred molecules, as 
is studied here. 

Correspondingly, MC simulations of the kind reported 
here can also easily be run for the number of MC steps re- 
quired to access equilibrium liquid properties. To demon- 
strate this, we show in Fig. [2] the P(p) histograms gen- 
erated in our simulations as a function of the computing 
time invested, for T = 245 K at a fixed value of the ac- 
tivity z*. We find that P(p) converges to well-defined 
values for all densities considered here. In addition we 
also check that the energy autocorrelation function de- 
cays to zero, in all SUS windows, in a number of steps 
smaller than the simulation length. 

We also note that Ref. [30] demonstrates that the 
dynamical properties of the liquid state of ST2 in the 
low density region near the optimal network density 
(0.83 g/cm 3 ) have a very simple and well-defined depen- 
dence on T, including in the region where T approaches 
T c and the LDL emerges as a distinct phase. 

V. ANALYSIS OF THE STRUCTURE OF THE 
LIQUID IN RELATION TO THE CRYSTAL 

Having shown in the previous section that liquid equi- 
librium is established in our simulations, we next test if 
crystal nucleation is occurring on the MC time scale of 
our runs. This is particularly relevant in view of the con- 
clusions presented in Ref. [23], in which it is argued that 
the LDL phase does not exist and that the only well- 
defined state of the system at low density is the crystal. 
We therefore must quantify the degree of crystalline or- 
der in our system at low T, especially in the LDL region, 
where crystallization might occur on a shorter time scale 
due to the similarity between ice and the liquid in terms 
of density and local structure. 

The degree of crystalline order can be quantified us- 
ing the Steinhardt bond order parameters [32] based on 
spherical harmonics of order I = 3 and I — 6, which 
are particularly suited for discriminating disordered fluid 
configurations from the open structure of hexagonal (as 
well as cubic) ice. For each particle we define the complex 
vector, 

Hm{i) = Jj-r- ^ Y lm(hj), (1) 

i=i 

where the sum is over the Nb(i) neighbors of particle i. 
Yi m is a spherical harmonic of order I and m, and r\y is a 
unit vector pointing from the oxygen atom on molecule i 
to that on molecule j. In the case of I = 3 two molecules 
are considered neighbors if their oxygen-oxygen distance 
is smaller than 0.34 nm, the position of the first minimum 
of the radial distribution function. In the case of I = 6, 
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FIG. 1: Consistency of results with previously published MD 
simulation data [251 130j . Shown are isotherms of the (a) po- 
tential energy and (b) pressure. Open symbols indicate data 
obtained from the individual SUS windows. Curves show pre- 
vious MD results, starting from T = 230 K and spaced every 
5 K. Curves are thicker for T matching the SUS simulations. 

to be consistent with Ref. [23], we assume Nb(i) — 4 and 
neighbors are defined as the four closest particles. The 
dot product, 

i 

c l = ( 2 ) 

771 — — I 

where, 

( 1 2 V /2 

qim(i) = qim{i)/ I Y \qi m (i)\ 2 ] (3) 

\m=-i / 

and q_i m {i) is its complex conjugate, determines the de- 
gree of orientational correlation between neighboring par- 
ticles i and j. Fig. [5] shows the distributions of cj 3 for 
I = 3 and I = 6, for several densities at T = 240 K. The 
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FIG. 2: Example of the convergence of the density of states 
during SUS MC simulations. Shown is P{p) for T = 245 K 
and z* = 1.46 x 10~ 4 ton" 3 at various times over 19 CPU-days 
of simulation time. 

distributions have a single broad peak at higher p, and 
become more bimodal at lower p with a peak forming 
near — —1. However, the distribution goes to zero at 
C3 = —1. The bimodal shape of the c 3 J distribution at 
low p is characteristic of liquids with well formed tetra- 
hedral networks [33]. The figure also shows the same 
distribution evaluated in a hexagonal ice configuration 
at the same T. The crystal is characterized by a large 
peak centered at ~ —1 and a smaller peak located 
approximately at c£ ~ —0.1, with 1/3 the area of the 
large peak [34] . In cubic ice, only the peak at ~ — 1 
is present. The distribution of Cg is rather density in- 
dependent, but again very different from the crystalline 
one. 

Ref. [23] specifically investigated the global metric, 

N 
i=l 

and the related m-independent rotational invariant, 

/ l \ 

= JV I S Q^nQL n ■ (5) 

\rn=-l ) 

Fig. [3^c) shows Qq at T = 240 K evaluated from the 
different windows of densities as well as for hexagonal 
ice. Again, the configurations sampled in our SUS sim- 
ulations do not show any sign of crystalline order. The 
value of Q 6 in the liquid state, expected to be zero in the 
thermodynamic limit, is always small (Q§ rs 0.06) and is 
one order of magnitude smaller than the crystal value. 

It is possible that a crystallite large enough to affect 
system properties, but too small to significantly affect a 
global measure of crystallinity such as Q§, may be present 
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FIG. 3: Characterization of crystallinity in our system at 
T = 240 K. Panel (a) shows the distribution of bond correla- 
tion dot products P(c^) obtained using order I = 3 spherical 
harmonics for different SUS windows spanning our density 
range. Legend labels indicate p in each window. Also shown 
is the distribution for hexagonal ice at p = 0.87 g cm -3 . Panel 
(b) is the equivalent of (a), only using I = 6. Panel (c) shows 
Qe as a function of p for our simulations, as well as for hexag- 
onal ice. Also shown is the average fraction of crystalline 
particles in our system, using criteria based on both I — 3 
and I = 6. 
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FIG. 4: (a) Probability distribution function of the density 
P(p) obtained from SUS MC simulations at activity values 
of 2* = 1.145 x 10~ 4 for T = 240 K, z* = 1.46 x 10~ 4 for 
T = 245 K, z* = 1.85 x 10~ 4 for T = 250 and z* = 2.88 x 
10~ 4 at 260 K (all in nm ). The corresponding values of 
P are respectively P = 223, P = 192, P = 172 and P = 
130 MPa. Below T = 250 K, bimodality of the distribution 
emerges, signaling the appearance of two liquid phases with 
distinct densities, (b) Isothermal compressibility as a function 
of density along isotherms, obtained using Eq.|6] The location 
of the peaks matches the data presented in Ref. [25] . 



in the system. To test for this, we estimate the number 
of crystal-like particles in the system n crys . We label 
particles as being crystal-like using two definitions. For 
I = 3, we define a particle as crystal- like if it has at least 
three dot products with neighboring particles satisfying 
Cg* < —0.87. This criterion has been used in nucleation 
studies of tetrahedral systems [331 ES] and in our system 
the value of -0.87 is rather generous, in that the distri- 
bution for the cubic and hexagonal crystals near c^f ~ — 1 
falls to zero before reaching -0.87. For I = 6, we similarly 
define a particle to be crystal-like if it has three neighbors 
with Cg J > 0.70, a low estimate for the value at which the 
crystal distribution crosses the liquid distributions near 
Cg rj 1. We plot the average fraction n CTys /N of crystal- 
like particles using I = 3 and I = 6 at T = 240 K in 
Fig. ^c). The largest value for n crys /iV ss 0.024 suggests 
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FIG. 5: Density distributions P(p) at low T obtained by a 
reweighting of the T = 240 K energy-density distribution 
P(p,E). The high density peak moves to higher p with de- 
creasing T, reflecting the behavior of a simple liquid, while 
the location of the low density peak stays fixed because the 
geometric constraints for a well-formed random tetrahedral 
network depend sensitively on the density. Note that the nu- 
merical noise (see e.g. the apparent minimum at p w 0.89 
g/cm 3 ) is progressively amplified reweighting at smaller and 
smaller T. 

that the average crystal cluster size, in our system of 
about 200 — 300 molecules, is not larger than roughly six 
molecules, again confirming the lack of crystallinity over 
the entire density range studied. Even if we do not re- 
sort to average values, the largest number of crystal-like 
particles that we ever observe in our system, considering 
all of the configurations we sample, is n crys = 19, found 
when using I = 3. 

Our observation that crystal-like particles are rare at 
all densities, including in the range of the LDL phase, 
demonstrates that the crystal nucleation process occurs 
on a much longer time scale than that required for liquid- 
state relaxation. Further, the fact that those small clus- 
ters of crystal-like particles that do occasionally form 
subsequently disappear, shows that a finite and non- 
trivial nucleation barrier separates the liquid phase from 
the crystal phase at all densities. This demonstrates that 
the liquid phases simulated here are associated with free 
energy basins that are distinct from that of the crystal 
phase. 

VI. DENSITY OF STATES 

Fig. |4Ja) shows P(p) for all studied temperatures, at 
the chemical potential for which the density fluctuations 
are maximal. The data show the onset of a bimodal 
distribution below T = 250 K, consistent with the exis- 
tence of a LL critical point, and suggests the occurrence 
of two thermodynamically distinct liquids phases associ- 



ated with well-separated free energy basins. The inset 
shows the same data on a semi-log scale, to highlight the 
ability of the SUS method to provide an accurate esti- 
mate of P{p) for over 50 orders of magnitude. 

From the energy-density probability density P{p, E) 
obtained from the SUS simulations, all possible thermo- 
dynamic quantities can be calculated, for all values of 
chemical potential (limited only by the noise quality of 
the data). It is interesting to determine the behavior of 
the isothermal compressibility Kt, which in terms of the 
fluctuation in the number of particles is, 

(iV 2 ) - (iV) 2 
k B T PN K T = V f {N) X , (6) 

where pn is the number density and ks is the Boltz- 
mann constant. Fig. |4^b) shows the compressibility at 
several T. As expected, and in agreement with previous 
estimates for the Kj- extrema locus based on MD sim- 
ulations |25| . the line of Kt maxima moves to smaller 
densities (and hence lower pressure) as T increases away 
from the critical point. 

The statistical quality of our estimate for P(p, E) al- 
lows us to predict P(p) down to approximately T = 225 
K via temperature reweighting. The reweighted P(p) in 
Fig. [5] shows two well-resolved peaks, with a very shal- 
low minimum between them, indicating that the free en- 
ergy barrier for the system to jump from the LDL to 
the HDL phase and vice versa, even in a system of only 
a few hundred particles, is becoming significantly larger 
than the thermal energy. In the inset to Fig. [5j we plot 
— In P(p) to highlight this growing free energy barrier. 
While the density of the HDL phase progressively in- 
creases on cooling, the density of the coexisting LDL 
remains essentially constant. This highlights that the 
thermodynamic stabilization of the LDL phase (i.e. the 
establishment of an equilibrium network of tetrahedrally 
bonded molecules) requires a very precise density to be 
achieved. This strong coupling between the density, local 
tetrahedral geometry, and free energy is at the heart of 
the physics of water. 

VII. DISCUSSION 

The results reported in this article are consistent with 
the previous calculations of Liu, et al. [22], and show 
that the evidence for the proposed critical phenomenon 
in ST2 water is independent of the way the long range in- 
teractions are modelled. In addition, our results are also 
consistent with the possibility that the free energy of this 
model, projected onto the density, at low T is character- 
ized by two basins, both of which correspond to disor- 
dered liquid phases. We also show that even though the 
location of the proposed LL critical point lies on the ex- 
tension of the liquid free energy surface that is metastable 
with respect to crystal formation, the time required for 
homogeneous nucleation is sufficiently long in our sys- 
tem to allow for local equilibration in phase space. In 
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sum, the evidence presented here continues to point to 
the existence of a LL phase transition in supercooled ST2 
water. 

To further test for behavior consistent with a LL phase 
transition, a number of additional questions should be 
addressed. In particular, in a finite-sized system, a bi- 
modal density of states can occur in a non-critical sys- 
tem if the correlation length of the order parameter (here, 
the density) exceeds the system size. A finite-size scaling 
analysis of the density of states, using a range of system 
sizes, would therefore provide an important test for the 
proposed LL transition in the ST2 model. It would also 
be valuable to test if the free energy barrier for crystal 
nucleation is strongly affected by the system size, and by 
the constrained cubic geometry of the simulation cell, in 
order to more fully understand the relationship between 
the time scale for equilibrating the liquid and for crystal 



formation. We recommend such studies for future work. 
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